# Purpose: Create regression table of missing minutes reporting (at plant level)
#          against predicted levels of emissions and treatment

reporting_data <- read_dta(paste0(EMISSIONS_DATA_IN, "/CEMS_2018-12-31 to 2021-04-03_StackDay_Balanced_HistCFs_Master337s_3sample2trunc.dta"))

all_merged <- reporting_data %>%
  mutate(tot_report_mins = if_else(is.na(tot_report_mins), 0, tot_report_mins)) %>%
  mutate(missing_mins = 24 * 60 - tot_report_mins) %>%
  filter(!D_interregnum) %>%
  group_by(gpcb_id) %>%
  # Average missing as portion of the entire day
  summarize(gpcb_id_avg_missing = mean(missing_mins)/(24*60)) %>%
  ungroup() %>%
  left_join(read_dta(paste0(EMISSIONS_DATA_OUT, "/potential_max_emissions.dta"))) %>%
  left_join(read_dta(paste0(BASELINE_DATA_OUT, "/BaselineCovariates.dta"))) %>%
  # Taken from table D1 in ETS paper, 8 TPH panel
  mutate(
    abate_fac = case_when(
      as.logical(D_esp) ~ (1 - .997),
      as.logical(D_scr) ~ (1 - .94),
      as.logical(D_bf) ~ (1 - .99),
      as.logical(D_cyc) ~ (1 - .8)
    ),
    abate_cost = case_when(
      as.logical(D_esp) ~ 69.09,
      as.logical(D_scr) ~ 10.99,
      as.logical(D_bf) ~ 8.00,
      as.logical(D_cyc) ~ 5.60
    ),
    pred_emissions = abate_fac * AverageUncontrolledMass12
  ) %>%
  janitor::clean_names()

options(modelsummary_format_numeric_latex = "plain")
options(modelsummary_factory_latex = 'kableExtra')
fixest::feols(gpcb_id_avg_missing~abate_cost*d_treatment, data=all_merged, vcov = "hetero") %>%
  list() %>%
  setNames( linebreak("\\emph{Dependent Variable:}\nShare of Day Not-Reporting", align = "c")) %>%
  modelsummary::modelsummary(stars=c('*' = .1, '**' = .05, '***' = .01),
                             estimate = "{estimate}{stars}",
                             coef_map = c(
                               "d_treatment" = "Treatment (=1)",
                               "abate_cost" = "Predicted Abatement Cost",
                               "abate_cost:d_treatment" = "Predicted Abatement Cost $\\times$ Treatment"),
                             gof_map = list(
                               list(raw="r.squared",clean="R$^2$", fmt=2),
                               list(raw="nobs",clean="Observations", fmt=0)),
                             output = "latex",
                             booktabs=T,
                             escape = F
  ) %>%
  as.character() %>%
  # Extract tabular
  str_extract(".*(\\\\begin\\{tabular\\}\\s*([\\s\\S]*?)\\\\end\\{tabular\\}).*", group = 1) %>%
  write_file(., paste0(EMISSIONS_TABS, "/Table_F4.tex"))